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We present a theoretical study of spontaneous imbibition in a slit pore using a lattice-gas model 
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I. INTRODUCTION 

Spontaneous imbibition, namely the rise of a liquid in a capillary tube and more generally in a porous solid, is 
a ubiquitous phenomenon that has received a lot of attention over the years because of its crucial importance in 
many industrial processes (oil recovery, ink printing, textile dyeing, etc..) as well as in agriculture and biological 
sciences. In spite of the complexity of real porous media (in particular their disordered structure[T]), it is often 
observed that the distance of penetration of the liquid inside the solid increases asymptotically as the square root 
of time, as predicted by the classical Lucas- Washburn (LW) equation [2] that describes the fluid behavior in a single 
capillary. On the macroscopic scale, the \fi law results from the balance between a constant capillary driving force 
due to the pressure drop across the liquid/vapor meniscus and an increasing viscous drag. In recent years, in relation 
to the rapid development of nanofluidic devices, the applicability of these macroscopic concepts in strongly confined 
geometries has attracted much interest, on both the experimental^ and theoretical[4- 10J sides. The main conclusion 
of these studies is that the LW equation still works in capillaries with diameters of a few molecular sizes, the prediction 
being almost quantitative in the case of a simple Lennard- Jones fluid 0]. This is quite remarkable as the continuum 
hydrodynamic description of the fluid is expected to break down at very small length scales. 

In the case of complete wetting, molecular dynamics (MD) and lattice-Boltzmann simulations [3HE] have also reported 
the presence of microscopic films moving ahead of the main capillary front and following also a yf law, but with a 
different prefactor than that of the main front. Such precursor films are commonly observed in spreading droplet 
experiments and several theoretical models have been proposed to reproduce their diffusive behavior [11]. In the case 
of imbibition, wetting films can significantly affect the dynamics of the gas-liquid interface [12] and it is generally 
expected that they reduce the viscous drag, which of course is a crucial issue in microfluidics. Some aspects of the 
problem, however, are not as well understood and deserve more systematic investigations. In particular, one would 
like to better understand the experimental conditions required for the appearance of the films, the influence of the 
solid wettability on their behavior, and how they compete with the meniscus for the liquid coming from the reservoir. 

These are complex issues, due to the coupling of hydrodynamic and non-hydrodynamic modes in the vicinity of the 
contact line. Ideally, a theoretical study should require a full atomistic description of the system, but the very different 
length scales involved in the problem make it computationally very demanding. An important simplification consists 
in adopting a coarse-grained (or mesoscopic) lattice-gas description, as is done in the lattice-Boltzmann method which 
indeed appears as an efficient tool for modelling capillary filling 5j-[7]. In the present work, we further simplify the 
problem by treating the dynamical evolution of the fluid in configurational space only, the effect of the hydrodynamic 
modes being simply incorporated in effective parameters: only diffusionlike mechanisms are explicitly included. This 
minimal model, however, does account for liquid conservation which is the key ingredient at the origin of the constant 
slowing down of the imbibition front. A similar approach underlies the recent phase-field models of imbibition in 
disordered porous solids [T] fT3] IT4]. The study of such systems is the next step on our agenda [15]. and its theoretical 
treatment indeed requires the above simplifications. Here, we describe spontaneous imbibition in a slit-pore by using a 
dynamic (lattice) mean-field theory (DMFT) which is essentially the mean-field version of the Kawasaki spin exchange 
dynamics. This theory is consistent with the treatment of thermodynamics at the mean-field level and it has already 
proven useful to model the dynamics of adsorption/desorption phenomena in nanopores [16ITT8] . 

Although we shall mainly consider the case of primary imbibition, when a liquid invades a dry capillary and precursor 
films appear in the complete wetting regime (which is the situation considered in the recent simulation studies [3] [5]), 
we shall also briefly investigate the case where the liquid advances over pre-existing thin films (prewetted capillary) 
which is encountered in many actual situations, in particular in experiments with nanoporous solids [3] 119] . 
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II. MODEL AND THEORY 
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FIG. 1: (Color on line) Schematic representation of the slit pore geometry used in the DMFT calculations. The pore has a 
length L is the z direction and a width H in the x direction. Periodic boundary conditions are applied in the y direction. The 
liquid reservoir (in blue) is located at the left side of the pore. 



We consider the case of a slit pore schematically represented in Fig. [T] Before describing the dynamic mean-field 
approach, we first briefly recall the macroscopic description that leads to the LW \fi law. As pointed out in the 
introduction, the dynamics of spontaneous imbibition under steady conditions (and neglecting evaporation, gravity, 
and incrtial effects Ij) is governed by the balance between the capillary driving force due to the Laplace pressure 
across the liquid/gas meniscus and the viscous drag of the liquid. In a slit geometry the meniscus is cylindrical and 
the Laplace pressure is given by 

= 2 7 , 3 cosg 

H 1 ; 

where H is the width of the slit, ji g is the surface tension of the liquid-gas interface and 9 is the contact angle between 
the liquid and the solid walls (assuming that in the long-time limit one can replace the dynamic contact angle by 
its static value). For an incompressible fluid, this gives rise to a constant gradient pressure over the whole imbibed 
region, 

where h(t) is the position of the meniscus inside the pore at time t (neglecting the internal structure of the interface). 
For a parabolic fluid- velocity profile normal to the walls (Hagen-Poiseuille laminar flow) , the viscous drag implies that 
the average velocity of the front is proportional to AP, 

h(t) = --AP (3) 
V 

where rj is the fluid viscosity and k = H 2 /12 is the pore permeability. This leads to the differential equation 



which after integration gives 



hh = 7lg Hc 0S e 

or? 



h{t) - h{t = Q) = { llME^lf/^. (5 ) 
6rj 



The progression of the liquid inside the pore is thus described by a diffusive law with an effective diffusion coefficient 
that scales like the pore width. Since the viscous drag is smaller in larger pores, the liquid progresses faster, as indeed 
observed experimentally. 

Our presentation of DMFT closely follows that of Refs. |I6lHH] and we refer the reader to these recent papers for 
more details (see also Ref. [20 and the comprehensive review [2 1 ). This mean-field version of the Kawasaki spin 
exchange model was first applied to the study of phase separation and surface enrichment in solid binary alloys 22 . 
More generally, it is appropriate to describe systems with a diffusionlike behavior. The fluid inside the pore is modeled 
as a simple cubic lattice gas with Hamiltonian 

n = ~ w ff nin i ~ Ws f ™ ( 6 ) 

<ij> i, surface 
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where n, is the occupancy variable at site i (n, = 0, 1) and the first sum runs over all nearest-neighbor pairs associated 
with fluid sites while the second sum runs over the fluid sites within the layers adjacent to the two walls (the pore 
walls being aligned with the (100) planes of the cubic lattice); Wff and w s f denote the nearest-neighbor fluid-fluid 
and solid-fluid attractive interactions, respectively, and the wettability of the pore walls is thus controlled by the 
parameter a = w s f /w / / . The corresponding bulk lattice gas has a symmetric coexistence curve pi (T) = 1 — p g (T) 
with chemical potential p sa t — — 3u>//, independent of temperature T. 

This type of lattice-gas model has been extensively used in the literature to study adsorption phenomena and 
phase transitions on substrates and in confined geometry (see e.g. [23-26 ). Here, for simplicity, the effect of the 
solid is felt only in the layers immediately adjacent to the pore walls, as in the recent applications of DMFT to 
capillary condensation and evaporation|16- 18 . This is of course a simplification of the actual physical situation, 
where one typically has a long-ranged van der Waals interaction which decays as a power-law with the distance from 
the surface [27 . Note also that we do not include any effect of surface corrugation in our calculation, although they 
may play an important role at low temperature. These two effects can be introduced in the theoretical treatment 
at the price of increasing numerical complexity. In the present set-up, the system is homogeneous in the y direction 
parallel to the walls and is effectively two-dimensional. 

Within the master equation approach|22). one can write down an exact equation for the temporal evolution of the 
mean site occupancy variable 16- 18] [20l I2T] 

Pi{t) = {n i ) t = ^n i P{{n},t), (7) 
{„} 

where P({n},t) is the probability to find the occupancy configuration {n} at time t. Assuming that the transport 
is due to a hopping process between nearest-neighbor sites corresponding to a Kawasaki exchange dynamics, this 
evolution equation can be written as 

|r = -£'«(*> (8) 

j/i 

where the summation runs over nearest neighbors of site i and the flux Jij (t) from site i to site j is given by 

Jij{t) = (Jij{{n}))t = ( w ij{{ n }) n i( 1 - Tij) - Wji({n})nj(l - 7i,-)) t , (9) 

with Wij({n}) the transition probability for transitions from site i to site j for a configuration {n}. Through the 
Monte Carlo method with Metropolis transition probabilities, one could obtain an explicit realization of this dynamical 
process. In the DMFT, one instead uses a mean-field approximation and merely replaces the occupancy variables 
by their ensemble average (which amounts in particular to assuming that the occupancy variables are dynamically 
uncorrelated) . This yields 

Jij(t) = WijPi(l - pj) - WjiPj(l - pi) , (10) 
and the evolution equation for the site densities becomes 



dp 

with transition probabilities 
where 



gt - - Kjpi(i - Pj) - Ujipji 1 - Pi)] (ii) 

j/i 



w ij({p}) = woexpi-PEij) (12) 



E- - i °' Ej < Ei (13) 



and 



E i = -Wff J2pj-\ 
j/i L 



if i is a "bulk" site ...... 

if i is a "surface" site. ^ ' 



The parameter wo is an elementary jump rate which sets the time scale. We take the same Wq for the "bulk" sites 
and the "surface" sites adjacent to the pore walls, thereby implying that the diffusion coefficient for a single molecule 



4 



does not depend on the solid-fluid interaction, which is an additional approximation iJSj ■ The above dynamics could 
also be interpreted as resulting from a coarse-grained approximation of the full hydrodynamic treatment. In such an 
interpretation, the parameter w incorporates in an effective way some average information about the hydrodynamic 
modes and could for instance be taken as proportional to the ratio of the permeability to the viscosity k/tj, with 
an explicit dependence on the pore width H (see above). However, in what follows, we neglect such remains of the 
hydrodynamic treatment and consider wq as independent of H. (As a result, we do not expect to recover the proper 
trend for the front velocity as a function of the pore width and shall not study this aspect of the problem.) 



As discussed in Refs. |16l 121) . Eq. (Ill may be recast into the discrete version of a Cahn-Hilliard [2"9"] or dynamic 
density functional theory[5ni EI] equation. The characteristic feature of the DMFT is that the mobility coefficient 
has an explicit expression and is related to the local site densities via the Metropolis transition probabilities. It 
is also worth noting that a coarse-grained continuity equation with a current J proportional to the gradient of a 
chemical potential (i.e. a generalized Cahn-Hilliard equation with a Ginzburg- Landau type free energy) is at the 
core of the phase-field models that have been developed recently for studying spontaneous and forced-flow imbibition 
in disordered porous solids [TJ Q31 [14]. In such a continuum (mesoscopic) approach, the dry and wetted states are 
considered as two "phases" of the system ((f) = ±1 where <j> is the locally conserved field) and the influence of disorder 
on the capillary force and the viscous drag are included in a locally random, but concentration- independent, mobility. 
Hydrodynamic effects are thus hidden in this effective quantity. 

An important feature of the DMFT, which will allow us to make a connection between the properties of the solution 
of the evolution equation and the adsorption isotherms in the same system, is that the site densities obtained from 
the minimization of the static Hclmholtz free energy ^"[{/Oj}] in the canonical ensemble (or from the minimization of 



the grand-potential Vt — T — p, J2i Pi m the grand canonical ensemble) are steady-state solutions of Eq. (111. Indeed, 
T in the mean-field approximation is given by 

PF[{pi}] = ln P l + ~ P*> ~ Pi)1 ~ @ W ff Y PiPi ~ P Ws f Y P l ' ( 15 ) 
i <*J> i, surface 

and the equilibrium (or at least metastable) pi's that satisfy the set of coupled equations 

-e*-= exp \P(n-E i )] (16) 



are also solutions of Eq. ( 11 ) when the left-hand side is set to zero [16] . 

In practice, we have integrated Eq. (11 1 via Euler's method (i.e. pi(t + At) = p,i(t) — At^j/i Jij(t)) by using the 



dimensionless time step woAt = 0.2. This value was small enough to ensure the stability of the solution in all cases 
and it actually corresponds to a very small variation of the site densities (which leads to a serious technical problem 
and makes the numerical calculations rather time-consuming). For the presentation of the results, we have used a 
larger time scale, to = 10 5 wo<At. 



III. RESULTS AND DISCUSSION 



Most of the results presented in this section correspond to the case of a "dry" capillary, where the liquid in the 
reservoir is initially in equilibrium with the vapor in the pore at p, = p sa t ■ This is realized by fixing the density of the 
fluid at the saturated liquid density pi(T) in the first (left) layer in the z direction and at the saturated vapor density 
p g (T) in the rest of the pore. Then, at time t — 0, we let the system evolve according to Eq. (11) and the liquid 
enters the pore. This is the case considered in the recent MD simulations^ [51 [5] where the capillary walls are initially 
taken lyophobic and then switched to lyophilic at t = 0. As in these simulation studies, the pore in our set-up has 
its right end closed by a wall that prevents the fluid from escaping. An alternative set-up consists in starting from 
a "gas" configuration inside the pore with the adsorbed films along the walls already formed. We shall also briefly 
consider this second situation as it is as well encountered in experiments with nanocapillaries[3l 119] . 

The calculations are done at the reduced temperature T* = fe^T/io// = 1 (the bulk critical temperature is T* — 1.5) 
for which p g 0.07 and pi ~ 0.93. Unless explicitly stated, the pore width H is fixed at 30 lattice constants (i.e. 28 
fluid layers). For comparison, some calculations for other values of H have also been performed. 

For future reference, it is instructive to first examine the influence of the wettability parameter a on the adsorption 
isotherms obtained from the free-energy, Eq. (15), in the grand-canonical ensemble. These are shown in Fig. [2] as a 
function of the relative activity \/X sa t = exp[/3(/x — p sa t)]- As usual, these isotherms have been obtained by starting 
from a low relative activity and increasing A in small steps. This type of step-like isotherms is well documented in the 
literature (see e.g. [26] ) and our primary interest here is in the two-dimensional layers of different thicknesses that 
form successively as a increases. Note that some of them are actually metastable since the capillary condensation (not 
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FIG. 2: (Color on line) Grand-canonical adsorption isotherms for T* = 1 in a pore of width H = 30 for various values of the 
wettability parameter a. From right to left: a = 0.8, 0.9, 1, 1.2, 1.6, 2 (for a = 0.8, the part of the isotherm that extends beyond 
bulk saturation is not shown). The figure focuses on the range of X/\ sa t associated to the formation of two-dimensional layers. 
The jump in density associated with capillary condensation is not shown. 



shown in the figure) occurs before the layering transition (for instance, this is the case of the second layering transition 
for a — 1.2). This distinction, however, is irrelevant for the dynamic evolution during spontaneous imbibition as will 
be seen below. 

As a necessary ingredient to the classical description of imbibition, it is important to know the variation of the 
static contact angle 9 with a. The contact angle is computed from Young's equation using the interfacial tensions on 
an infinite planar surface obtained from the mean-field free energy in Eq. (15). We refer the reader to Ref . [32] for 
more details on the calculation and to Fig. 2 in Ref.fTS] for the results at T* — 1. (Note that cos# = for a = 0.5 
due to the symmetry of the nearest-neighbor lattice-gas Harniltonian[2H H3; note also that the weak dependence of 
6 on the direction of the interface [34] is neglected in this calculation.) At this temperature, the wetting transition as 
a function of a is first-order and occurs slightly below a = 0.9. If one neglects the effect of confinement, this explains 
the behavior of the a — 0.8 isotherm in Fig. [2j with no layering transition before A = \ sa t- 



A. Dry capillary 

We first discuss an important methodological issue that must be considered in numerical studies of an initially dry 
capillary. It concerns the length L of the system. As usual, L must be chosen large enough to avoid undesirable finite 
size or edge effects. When precursor films are present, it turns out that this is a much more severe requirement than 
for standard calculations of adsorption isotherms. This is illustrated in Fig. [3] where we compare the time evolution 
of the average fluid density profile p{z) for a — 2 and three different pore lengths, L = 2000, 4000 and 6000. The 
characteristic features of the profile will be discussed in more detail below but one can see at once that the advance of 
the main interface is modified as soon as the average fluid density at z = L becomes significant (typically p(L) > 0.1) 
due to the fact that the precursor film has already reached the pore end: this results in a significant increase in the 
velocity of the interface. The figure shows that for L as large as 4000 lattice constants (which is eventually the length 
used in our calculations), the numerical results are biased if the imbibition time exceeds t ss 400 to. Note that the 
main interface at this time has only reached a distance z pa 86, much smaller than L ! The length of the capillary is 
therefore an issue that must be treated very carefully in this type of calculation [35]. 

We now discuss the main results of our calculations. Fig. [4] shows the time evolution of the average fluid density 
profile p{z) in the capillary for different values of a. 

In all cases, the observed equidistance between the successive profiles at times t/t = 1, 2 2 , 3 2 , 4 2 , ...20 2 shows that 
the main front, which corresponds to the largest variation in the average density and therefore can be associated to 
the liquid-gas interface in the center of the pore, advances like \ft (see also Fig. 5). In the complete wetting regime 
(a = 1, 1.2 and 2), the main front is preceded by a smaller front that advances faster, but also with a \[t law (for 
a = 1.2 there are actually two small fronts, the first one advancing faster than the second, see below). These fronts 
can be associated to thin precursor films propagating along the pore walls. For a — 0.8, which corresponds to a partial 
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FIG. 3: (Color on line) Profiles of the average fluid density p(z) in a pore of width H = 30 as a function of time for a = 2 and 
three different system lengths L — 2000 (black solid lines), L — 4000 (red dashed lines) and L — 6000 (blue dot-dashed lines). 
From left to right, the profiles correspond to t/to = 4 2 , 8 2 , ...20 2 . The distance z is put on a logarithmic scale, so that z = 4000 
and z — 6000 are virtually indistinguishable on the scale of the figure. 
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FIG. 4: Profiles of the average fluid density p(z) in a pore of width H — 30 (L — 4000) as a function of time for different 
values of a: (a) 0.8, (b) 1, (c) 1.2, and (d) 2. The successive curves correspond to t/to = 1, 2 2 , 3 2 ...20 2 . 



wetting situation, there is no precursor. On the whole, these observations are in line with the results of the recent 
off-lattice MD and lattice-Boltzmann simulations 4.5. The dependence on a, however, is nontrivial and deserves a 
more detailed investigation. 

There is indeed an interesting difference between the evolution with time of the total fluid uptake and that of the 
position of the meniscus. This is illustrated in Fig. [5] for values of a ranging from 0.7 to 3. The uptake (per unit 
lateral surface) is defined as T(t) = J[p(z,t) — p(z,t = 0)] dz, where we have subtracted the initial number of fluid 
molecules as we are only interested in the filling process for t > 0. One can see in Figure [5^, that the uptake increases 
with a and appears to saturate for large values of a (well above the transition to complete wetting). The fit to a 
\Ji law is excellent so that one can define an effective diffusion coefficient Dr — T(t) 2 /(t/to) whose variation with a 
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FIG. 5: (Color on line) (a) Total fluid uptake in the capillary as a function of the square root of imbibition time for various 
values of a (from bottom to top: a = 0.7,0.8,0.9, 1, 1.2, 1.6,2,3. The solid lines correspond to a s/t fit. (b) Advance of the 
liquid-gas meniscus. The symbols are the same as in (a) but the order of the curves is different (at the scale of the figure the 
results for a — 0.8 (circles) and 1 (squares) are almost indistinguishable). In both figures, the transient behavior for t/to < 1 
is not shown. 



is shown in Fig. [6^i. On the other hand, the behavior displayed in Fig. [SJd for the advance of the meniscus appears 
at first sight surprising. Although the results can be again very well fitted to a \fi law, the variation of the velocity 
with a is not monotonic. The effective diffusion coefficient defined by D — h(t) 2 / (t/to) first increases (from a = 0.7 
to ps 0.88), then decreases, and finally saturates, as shown in Fig. [6b. 




).5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 

a a 



FIG. 6: (Color on line) Effective diffusion coefficients Dr and D for the time evolution of the total uptake (a) and the advance 
of the meniscus (b) as a function of the wettability parameter a. In (b) the red dashed line corresponds to the formula 
D — D(a — 0.88) cos 8 (with 6(0.88) w 0). In (c) the product (DDr) 1 ' 2 is shown to become approximately constant in the 
complete wetting regime, in agreement with Eq. (17 1. 



It is clear that the first regime in which the velocity of the meniscus increases with a coincides with the situation 
of partial wetting. Moreover, for 0.7 < a < 0.88, the effective diffusion coefficient is reasonably proportional to 
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cos6, in agreement with the macroscopic picture for the driving capillary force. For smaller values of a, however, 
the relation between D and cos 9 is no more linear, as if the actual driving force were smaller; since we only consider 
here the asymptotic steady state and no variation of the contact angle with velocity is observed (apart from small 
fluctuations |36j ) . the concept of dynamic contact angle cannot be invoked to explain the observed deviations [3 7j. They 
are more likely related to the deformation of the meniscus which can be seen in Fig. [7] where some typical interface 
profiles at the same imbibition time are shown. The interface is here defined as the density isocontour p(x, z) = 0.5 
(the intrinsic width of the interface being equal to a few lattice constants, which is standard at this temperature). For 
a > 0.7, the menisci can be well fitted to a semi-circle with radius R and the contact angle extracted from the tangent 
at the wall (i.e. cos (9 « H/2R) is in reasonable agreement with the static angle computed on a flat interface [TH]: see 
the caption of Fig. 7. On the other hand, for a < 0.7, some flat portion appears in the center of the pore, a feature 
which may possibly be attributed to lattice artifacts that become more and more significant as 9 approaches ir/2. 




FIG. 7: Meniscus profiles for different values of a at t/to = 250. From bottom to top a — 0.7, 0.8, 1, 1.2. For clarity, the curves 
are arbitrarily shifted vertically by a constant and the wetting layers adjacent to the walls (when present) are not shown. The 
solid lines correspond to a semi-circular fit. From this fit, we obtain cos# = 0.48, 0.83 for a = 0.7, 0.8, to be compared with the 
values, 0.54,0.805, computed on a flat interface [18]: in both calculations, cos# ~ 1 for a > 0.9. 



The unexpected decrease of D with a in the complete wetting regime is related to the presence of precursor films [42 j. 
Indeed, because of mass conservation, the velocities of the meniscus and the precursor front are not independent. A 
semi-quantitative argument is developed in the Appendix where we investigate the coarse-grained behavior of the 
DMFT evolution equation on a scale where the structure of the interface may be neglected. In particular, we derive 
the relation 



y/DD r oc Sfj, (17) 

where Sfj, = fi sa t — fi* and fj* is the value of the chemical potential at the main interface. 5 ft is a small quantity 
that does not vary much in the complete wetting regime (Sfi/wft ~ 0.023 ± 0.002) and is approximately equal to 
(pi — Pg)jiv/R, in agreement with the macroscopic Gibbs-Thomson equation. As shown in Fig. pp, the product 
y/DDr is indeed constant in the whole range a = 1 — 3 with an error less than 5% (for a = 0.9 the product is smaller, 
possibly because the complete wetting regime is not fully established). It is remarkable that the product of the two 
effective diffusion constants D and Dp is only controlled by Sfj,. Since the speed of the whole imbibed fluid, including 



the precursor film and beyond, increases with a (which is expected on physical ground), Eq. (17 1 tells us that the 
speed of the meniscus must decrease. 

We now focus on the dynamics of the precursor films. Figure [8] shows the position of the precursor films as a 
function of y/t/to for increasing values of a. As was observed in Fig. [4] there may exist one or two moving fronts: 
h\(t) (resp. h,2(t)) denotes the position of the front corresponding to the filling of the first (resp. second layer) 
adjacent to the pore walls (the fronts are defined after averaging over the transverse direction x and thus incorporate 
the layers adjacent to the two walls and the gas in between). The first front is not well defined for a > 1.3. One can 
see that the velocity of the fronts increases with a and that the first one develops a faster dynamics when both fronts 
coexist. (The motion is much faster than that of the meniscus, which is why we were forced to consider very large 
systems.) In all cases, the fit to a \fi law is excellent, in agreement with the coarse-grained analysis performed in the 



Appendix (see Eq. (A16|). 

It is important to note that the precursor films do not advance on completely dry walls. Indeed, as soon as one 
lets the system evolve, there is a very quick redistribution of the gas inside the pore due to the attraction exerted 
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FIG. 8: The positions h\(t) and /i2(t) of the precursor films as a function of the square root of time for several values of a in 
the complete wetting regime. For hi(t): a — 0.9, 1, 1.1, 1.2, 1.3 from bottom to top; for h^: a = 1.2, 1.3, 1.6, 2 from bottom to 
top. The solid lines correspond to a \ft fit. 



by the walls. This transient process has no relation with the progression of the liquid from the reservoir and it thus 
establishes a new inhomogeneous "canonical" equilibrium (accordingly, the chemical potential at z = L becomes lower 
than p sat as shown in Fig. 13 in the Appendix). This effect is not visible in Fig. [1] which only displays the average 
fluid density p(z), but it can be seen in Fig. [9] where we plot the time evolution of the density pi(z) in the first layer 
for a = 1 and 1.2. For t/t = 10, one finds /^(lOOO) « 0.25 for a — 1 and « 0.36 for a = 1.2, which are both values 




(a) 




2(100 1000 2000 3000 



FIG. 9: Density pi(z) in the first layer for a = 1 (a) and a = 1.2 (b) at times t/to = 10,40,90, 160,250,360. In (a) the large 
and small jumps are associated to the passage of the precursor front and the passage of the meniscus, respectively. In (b) there 
is an intermediate jump associated to the passage of the second front displayed in Fig. rj;. 



significantly larger than the density p g ss 0.07 that was imposed at the beginning of the calculation. As expected, 
this initial density increases with a and we have checked that it also increases with the pore width as there are more 
particles than can be "pumped" from the center of the pore to be adsorbed on the walls. This inhomogeneous "gas" 
configuration should thus be considered as the actual initial state of the pore before imbibition begins. 

Fig. [9] shows that the filling dynamics of the first layer can be divided into distinct stages. For instance, for a = 1, 
there is first a gradual filling, then a large jump corresponding to the passage of the precursor front, a new gradual 
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filling, and finally a small jump corresponding to the passage of the meniscus. (For a = 1.2 the first stage is absent but 
there is an additional intermediate jump that corresponds to the passage of the second front ahead of the meniscus 
seen in Fig. [3J; .) The first stage of this process is analyzed in Fig. 10 which shows the time evolution of pi(z) at a 
fixed distance z for a = 1. The gradual filling of the layer is well fitted by the formula pi(z, t) = a + by/t/z 2 so that 
this process is also diffusive. 




FIG. 10: Time evolution of the density in the first layer at a fixed distance from the pore entrance (from left to right, 
z = 400, 800, 1200, 1600) for a = 1. The solid lines correspond to the fit p x (z, t) « 0.221 + 8.23y/{t/t )/z 2 . 



Remarkably, the values of the average density p in the pore before and after the passage of the precursor front are 
in good agreement with the values of the density observed before and after the corresponding layering transitions in 
the adsorption isotherms : for example, these densities in Fig. [4] are respectively 0.083 and 0.110 for a = 1, to be 
compared to 0.087 and 0.115 in Fig. [5] for a = 2, the density jumps from 0.148 and 0.180 at the passage of the front to 
be compared to 0.152 and 0.185 before and after the second layering transition in Fig. [2j The corresponding chemical 



potentials are also very close (for instance, pi/wff f=s —3.11 in Fig. 13d and pi/wff = —3 + T* ln(Ai/A sa t) ps —3.10 
in Fig. [2] for a = 1). This illustrates the fact that the same states are encountered during imbibition and adsorption 
processes, the possible metastability of a state with respect to capillary condensation being irrelevant dynamically. 

One could thus interpret the jump associated to the passage of the precursor front as a dynamic layering transition 
that occurs when p (or pi) reaches a certain critical value (for example, p\ ss 0.34 for a = 1). In this interpretation, 
the speed of the precursor is directly related to the time needed to reach this critical density. Since the initial density 
in the layer increases with the pore width, as noted above, this time gets shorter and the precursor moves faster as 
H increases, which is indeed observed in our calculations. Moreover, the intermediate stage corresponding to the 
gradual filling of the layer may be totally absent. Knowing the initial configuration of the fluid in the pore is therefore 
important to understand the dynamical behavior of the film, which is itself related to the dynamical behavior of the 
meniscus, as shown by Eqs. (A17) in the Appendix. (Of course, the viscous drag also varies with the pore width, 
an effect which is not included, even in an effective manner, in the present treatment where we have taken wq as 
independent of H.) 



B. Prewetted capillary 

We now consider the situation of a prewetted capillary. The pre-existing films coating the pore walls were created 
by gradually increasing the relative activity X/X sa t from a low value up to 1, and one or two monolayers were thus 
obtained depending on the value of a, as shown in Fig. [2j Of course, these thin films are metastable since the true 
equilibrium situation at p sat corresponds to the pore entirely filled with liquid and, in any case, capillary condensation 
occurs before p sa t- This feature, however, is not important for the present purpose and this situation may possibly 
occur in actual experiments with nanoporous disordered solids, depending on the timescales involved in capillary 
condensation and imbibition[19 . 



Figure 11 shows the time evolution of the average fluid density profile p(z) for various values of a. Comparing with 
Fig. |4j one observes the absence of precursor fronts, as expected. The advance of the interface again follows a \fi 
law, but the speed of the meniscus now increases with a like the total fluid uptake, as shown in Fig. [12} Moreover, 
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both quantities saturate more rapidly with a than in the presence of precursor films. The whole dynamical behavior 
is therefore quite simpler. 
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FIG. 11: Same as Fig. |4]for the prewetted capillary. 




FIG. 12: Same as Fig. [6] for the prewetted capillary. 



It is generally expected that the meniscus moves faster in a prewetted capillary as the frictional force is reduced [15], 
Comparing Figs. [IT] and 12 with [5] and [6] shows that this behavior is also observed in our calculations, but this is now 
due to the fact that the meniscus no more competes with the precursor film for the liquid coming from the reservoir. 
The diffusion coefficient for the total uptake on the other hand is sm aller than in the previous set-up (note that the 
product y/DDr is almost the same in the two cases because Eq. (17) still holds and Sfi does not really change). 
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IV. SUMMARY AND CONCLUSION 



To summarize, we have investigated spontaneous imbibition in a narrow slit pore in the framework of a dynamic 
lattice mean-field theory that explicitly includes liquid conservation but neglects momentum transport. The theory is 
able to reproduce the main characteristic feature of imbibition in the quasistatic regime, namely the constant slowing 
down of the interface described by the macroscopic \ft Lucas- Washburn law. When the capillary is initially dry and 
the solid-fluid interaction is sufficiently attractive (so that the static contact angle is zero), we have found that thin 
precursor films move ahead of the meniscus, following also a \[t law, as also reported by recent simulation studies [3J- 
|B]. The jump in the average fluid density in the pore that is associated to the passage of the precursor front may 
be interpreted as a dynamic layering transition, whose existence can be predicted from the corresponding adsorption 
isotherms. Due to mass conservation, the meniscus and the precursor compete for the liquid coming from the reservoir, 
which induces a nontrivial dependence of the speed of the interface on the strength of the solid-fluid interaction. On 
the other hand, the total fluid uptake increases monotonously. This contrasts with the simpler behavior observed in 
a capillary where thin films adsorbed on the pore walls are pre-existing: in this case, both the speed of the meniscus 
and the fluid uptake increase with the solid-fluid interaction and rapidly saturate in the complete wetting regime. 

The extent to which the above predictions may be tested in actual experiments or in simulations remains an open 
question. We have emphasized that the system length is an issue that must be carefully treated in simulations when 
precursor films are present. At a different level, since our treatment has drastically simplified the complexity of the 
problem, the phenomena revealed by our calculations may be hidden or significantly altered. For instance, we have 
neglected a possible dependence of the molecular friction coefficient on the strength of the solid-fluid interaction in 
the first layers adjacent to the pore walls. This dependence has been invoked to explain the nonmonotonic variation 
of the film spreading rate on solid substrates with different surface energies [28]. This could be included in the present 
theory by modifying the jump rate parameter wq in the first layer (for instance, by taking w\ = wq exp(— Ca) where 
C is some constant). By reducing the speed of the precursor for large a, this would in turn accelerate the speed of 
the meniscus. The absence of hydrodynamics is of course another very serious limitation if one wishes to go beyond 
a mere account through an effective elementary time scale. In principle, momentum transport could be incorporated 
at a coarse-grained level by coupling the density to a velocity field satisfying the appropriate Navier-Stokes equation 
(see e.g. [3H1I1S])- This, however, would considerably complicate the numerical treatment and make the study of 
more complex geometries encountered in disordered porous media virtually impossible. 
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Appendix A: Coarse-grained description of the DMFT evolution equation for imbibition 



We start with the DMFT evolution equation, Eq. (11), namely, 

-gjr = -^2 Kj'M 1 _ Pi) ~ uaPii 1 - Pi)] ■ ( A1 ) 

i/i 

Due to the symmetry, the y direction is irrelevant in the problem which is effectively two-dimensional in the (x, z) 
plane. We consider the sum of the /3,'s over the x direction at z constant. It is then easy to show that for all nearest 
neighbors of i in the z = constant plane, 

^ [wijpi(l - Pj) - WjiPj(l - p^] = . (A2) 
<ij>, z constant 

Then 

d r 

Ot^PiX'Z) = - ^2{[ w x(z,z + 1)(1 - p(x,z + 1)) +w x (z,z- 1)(1 - p(x,z- l))]p(x,z) 

X X 

- [w x (z + 1, z)p(x, z + 1) + w x (z - 1, z)p(x, z - 1)](1 - p(x, z))\ (A3) 
where we have defined w x {z 1 z') = Wij with i — (x, z) and j ~ (x, z'). 
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On physical ground, we expect that the energy E(x, z) increases and the density p(x, z) decreases with increasing 
z at constant x. Considering the definition of the rates, it is then found that 

w x (z, z - 1) = w , w x (z, Z + l)= WQe -Pmx,z+\)-E(x,z)] 

W X (Z + 1,Z)= W , W X (Z -1,Z) = Woe -PlE(.x,z)-E(x,z-l)] K A V 

We define the local chemical potential p,(x, z) as an extension of Eq. (16): 




FIG. 13: Chemical potential profile in the center of the pore (solid line) and in the first two layers adjacent to the pore walls 
(dashed and dotted-dashed lines) for t/to = 250 and a = 0.8 (a), 1 (b), and 1.2 (c). On the scale of the figure the three curves 
are indistinguishable. 



Mx,z)=ln P } X f +pE(x,z) . (A5) 
1 p( x, z) 

Then, 

| P(x, z)=-w J2 {(e-^*+V-^ - l) P (x, z + l)(l- p(x, z)) 

X X 

_ ^-^(x,z)-„(x,z-l)] _ _ p ^ z _ i))| . (A6) 

Except possibly in the vicinity of the interfaces, the chemical potential is expected to change smoothly and slowly. 
To take the continuum limit which corresponds to a phase-field model, one has to reinstall the lattice spacing "a" and 
take the appropriate limit a — > 0. In any case, (3[p,(x, z + a) — p(x, z)] <~ a <C 1, which leads to 

d ( 

-Q t X! P( x > z ) = P w o ^2 { z + a )- z)]p{x, z + a)[l- p{x, z)] 

X X 

- \p(x, z) - p(x, z - a)]p(x, z) [1 - p(x, z - a)} j 

= /3w X] {{K x > z + a ) + K x jZ~ a ) - 2 M#) z)]p{x,z)[l - P( x i z )} 

X 

+ \fj,{x, z + a) - p,(x, z)](l — p(x, z))[p(x, z + a) - p(x, z)] 

- [p(x, z) - p(x, z - a)]p(x, z)[p(x, z) - p(x, z - a)] \ . (A7) 
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In the continuum limit, a — > and wqcl 2 — > wq, one finds 

P( x > z ) = X! {^(^ ^)p(a;, z) [1 - p(x, z)] + d z p(x, z)d z p(x, z) [1 - 2p(x, z)] } 

X X 

= Pw ^ \d 2 z p{x, z)p(x, z)(l - p(x, z)) + d z p(x, z)d z [p(x, z)[l - p(x, z)]] } 

X 

= f3w <9 Z { ^2[d z p(x,z)]p(x 1 z)[l - p(x,z)}\ (A8) 

X 

When assuming that p(x, z) is independent of x, the evolution equation for the average density profile p{z) = 
(1/H) J2 X p( x i z ) can finally be written as 

^p{z) = d z (M{z)d z p{z)) (A9) 

where M(z) = Pwq{\/H) J2 x p( x : — p( x > z )} 1S a local, density-dependent mobility (the time dependence of p(z) 
and p(z) is left here implicit). A phase-field description is thus recovered at a coarse-grained level because the variation 
of the local chemical potential must be smooth enough to allow the expansion of the exponential and more generally, 
the variation of all z-dependent quantities must be smooth enough to allow the expansion in the lattice spacing a. 
Note that here p{z) is a true fluid density profile and not simply a mathematical object identifying the phases as in 
usual phase-field models. 

We now consider the physical situation encountered in the present study, in which the rapid variations of p and 
p are localized at h(t) (the position of the main front) and, when precursor films are present, at hi{t) > h(t). For 
simplicity, we assume that there is only one precursor front as found for a = 1 and a — 2; the calculation can be 
easily generalized to the situation where two or more fronts are present. We use the above phase- field description, 
neglecting the details of the fluid near the interfaces that are assumed to be sharp, and we focus on the region behind 
hi(t). For h(t) < z < hi{t) the averaged density p{z) = (l/H) J2 X p( x , z ) thus accounts for the two films on the pore 
walls and the gas in between. Inspection of the numerical solution shows that (i) p(x, z) is indeed almost independent 
of x far from the interfaces, as can be see in Fig. |13[ (ii) the mobility M(z) is approximately piecewise constant in this 
region, M(z) rss Mi = (3w pip Q for < z < h(t) and M(z) M* > Mi for h(t) < z < hi(t). Reinstalling explicitly 



the time dependence in Eq. (|A9|), we thus have 

i) 



Mi JW«,t) - < z < h(t) 



dt 



p(z,t) = { 2 (A10) 

M* &rn(z,t) , h(t) <z< hi(t) . 



At zeroth-order, a solution of this equation is provided by a piecewise description with p(z,t) linear in z and p(z,t) 
constant, 

= ^ sat ~ {tisat ~ ) for < z < h(t) (All) 

P(z,t)=pi J 

and 

p(z, t) = p* - (p* - Pi) hl z it} h X{ t) | for h (Q < z < hl (Q (A12 ) 
p{z,t)=p* i 

where p* , p* and p\ are considered as empirically fitted parameters (they can actually be evaluated by considering the 
adsorption isotherms, see Fig. 2). Note that this description is of course approximate and that corrections are present. 
For instance, as seen by blowing up the density profiles in Fig. 4, one finds that the density very slightly decreases 
between and h(t) (as an immediate consequence of the decrease in p) as well as between h(t) and hi(t). There is 
also a linear variation of the mobility M(z). These corrections however are small and well within the uncertainty of 
the present coarse-grained picture. The situation is more serious in the region ahead of the most advanced front (h±(t) 
if there is a precursor film, h(t) otherwise). Due to the attraction exerted by the pore walls, the chemical potential 
and the density are not simply those of the homogeneous gas and, as illustrated in Fig. [13J the chemical potential 
significantly varies with z and strongly deviates from a linear behavior (as is also observed for the density and the 
mobility). Computing p{z) and p(z) in this region is out of the scope of the coarse-grained description but one can 
reasonably assume (and check numerically) that these quantities scale like h(t) or hi(t), for instance p{z) = p(z / hi(t)) 
when there is a precursor film. 
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To determine the time evolution of the density profile for < z < h\(t) we then integrate Eq. (A10) over z in 
infinitesimal intervals around h(t) and h\(t). This yields 



h(t)+e 



h(t)-e 



dp(z,t) 
dt 



dz = M 



Mi 



, dp(z, t) , 
dz 

Psat — P 



h{t)+e 



Mi 



dp,(z, t) 



M* - 



Oz 

p* - fj,x 



\h(t)-e 



h{t) h x {t)-h{t) 
and thus, using the Heaviside form of the density profile at z = h{t), 

(pi — p*)h(t) — Mi ^ sa * ~ ^* — M* M *~ Ml 



h(t) 



Similarly, integrating between hi(t) — e and h\{t) + e: 

(P* - P**) 'h\{t) = M* 



P - Pi 
hi(t) - h(t) 



hi(t) - h(t) 



M* 



hi(t) ' 



(A13) 



(A14) 



(A15) 



where p** , M**, and pf{l) are the average density, the mobility, and the derivative of p(z/hi) at z = h\{t) + e, 
respectively. The solution of these two coupled differential equations is easily found to be of the diffusive form, 



(A16) 



h(t) = VDt 
h x {t) = t/DJ , 



with D and D\ solutions of the equations 

{ Pl - P*)Vd = 2M l 
(p* - p**)a/dT = 2M* 



Psat - P 

, p* - pi 



2M* 



P - Pi 



£>i - VD 

«A'(l) 



_ _ ,-2M _ 
D 1 -VD y/Dl 



(A17) 



where £>i > £>. Once the various parameters are known, D and Z?i are easily computed. 

Finally, the time evolution of the total fluid uptake T(t) = L L dz [p{z) — p(z. t — 0)] is obtained from 



r(f) 



dp(z,t) , n/ rf \ d ^ z \ nf Psat-p* 
dz — - — = ~[M{z)— — J z=0 = Mi- 



di 



h(t) 



(A18) 



where we have used the fact that the derivative of p(z) vanishes at the right end of the pore when L is large enough. 
This leads to 



T(i) = y/Dri 



with 



'D T = 2Mi 



Psat — P 



D 



(A19) 



(A20) 



Eqs. (A16), (A19) and (A17), (A20) provide the full solution of the (simplified) problem. 
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